# Clear memory
rm(list = ls())

# Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  knitr, kableExtra, data.table, tidyverse, tidylog
)
options("tidylog.display" = NULL)
transf <- function(x) (exp(x - 1) - 1) * 100

##########################
##### BASE WELFARE TABLE
##########################

cross <- fread("data/simulation_fips_crosswalk.csv")

# welfare results: 1997
welfare_base <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:unemp_amenity, ~ transf(.x))) %>%
  as_tibble()

# welfare results: no geography
welfare_pe <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-off_mobility-off_mktacc-off_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:unemp_amenity, ~ transf(.x))) %>%
  as_tibble()

# welfare results: first-best
welfare_first <- fread("data/counterfactual/welfare_total_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  bind_rows(
    fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv")
  ) %>%
  mutate(across(total:unemp_amenity, ~ 1 / .x * lead(.x))) %>%
  mutate(across(total:unemp_amenity, ~ 1 * transf(.x))) %>%
  as_tibble() %>%
  slice(1)

nonattainment <- fread("data/greenbook/nonattainment_status.csv") %>%
  filter(year == 1997) %>%
  filter(fips %in% cross$fips) %>%
  distinct(fips, nonattainment)

workers <- fread("data/payroll_cf.csv") %>% as_tibble()

# total number of workers in billions
workers_total <- sum(workers$employment, na.rm = T) / 1e9

welfare_base_industry <- fread("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  select(fips, industry, welfare, cf_labor, cf_wage, real_wage, cf_c_price, amenities) %>%
  # consumption equivalent change in welfare (100%) * base consumption * (1%/100%) * (10000 $ / 1$)
  mutate(
    total_wage_gain = transf(welfare) * cf_wage / cf_c_price / 100 * 1e4,
    cons_wage_gain = transf(real_wage) * cf_wage / cf_c_price / 100 * 1e4,
    amenity_wage_gain = transf(amenities) * cf_wage / cf_c_price / 100 * 1e4,
  ) %>%
  group_by(industry) %>%
  summarise(
    total_income_gain = sum(cf_labor * total_wage_gain * workers_total),
    cons_income_gain = sum(cf_labor * cons_wage_gain * workers_total),
    amenity_income_gain = sum(cf_labor * amenity_wage_gain * workers_total),
    wages = sum(cf_wage / cf_c_price * 1e4 * workers_total * cf_labor, na.rm = T)
  ) |>
  mutate(total_income_gain = cons_income_gain + amenity_income_gain)

welfare_first_industry <- fread("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  select(fips, industry, welfare, cf_labor, cf_wage, real_wage, cf_c_price, amenities) %>%
  # consumption equivalent change in welfare (100%) * base consumption * (1%/100%) * (10000 $ / 1$)
  mutate(
    total_wage_gain = transf(welfare) * cf_wage / cf_c_price / 100 * 1e4,
    cons_wage_gain = transf(real_wage) * cf_wage / cf_c_price / 100 * 1e4,
    amenity_wage_gain = transf(amenities) * cf_wage / cf_c_price / 100 * 1e4
  ) %>%
  group_by(industry) %>%
  summarise(
    total_income_gain = sum(cf_labor * total_wage_gain * workers_total),
    cons_income_gain = sum(cf_labor * cons_wage_gain * workers_total),
    amenity_income_gain = sum(cf_labor * amenity_wage_gain * workers_total), 
    wages = sum(cf_wage / cf_c_price * 1e4 * workers_total * cf_labor, na.rm = T)
  ) |>
  mutate(total_income_gain = cons_income_gain + amenity_income_gain)

welfare_pe_industry <- fread("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-off_mobility-off_mktacc-off_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  select(fips, industry, welfare, cf_labor, cf_wage, real_wage, cf_c_price, amenities) %>%
  # consumption equivalent change in welfare (100%) * base consumption * (1%/100%) * (10000 $ / 1$)
  mutate(
    total_wage_gain = transf(welfare) * cf_wage / cf_c_price / 100 * 1e4,
    cons_wage_gain = transf(real_wage) * cf_wage / cf_c_price / 100 * 1e4,
    amenity_wage_gain = transf(amenities) * cf_wage / cf_c_price / 100 * 1e4
  ) %>%
  group_by(industry) %>%
  summarise(
    total_income_gain = sum(cf_labor * total_wage_gain * workers_total),
    cons_income_gain = sum(cf_labor * cons_wage_gain * workers_total),
    amenity_income_gain = sum(cf_labor * amenity_wage_gain * workers_total),
    wages = sum(cf_wage / cf_c_price * 1e4 * workers_total * cf_labor, na.rm = T)
  ) |>
  mutate(total_income_gain = cons_income_gain + amenity_income_gain)

welfare_base_nonattain <- fread("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  left_join(nonattainment) %>%
  select(fips, welfare, cf_labor, cf_wage, real_wage, amenities, cf_c_price, nonattainment) %>%
  mutate(
    total_wage_gain = transf(welfare) * cf_wage / cf_c_price / 100 * 1e4,
    cons_wage_gain = transf(real_wage) * cf_wage / cf_c_price / 100 * 1e4,
    amenity_wage_gain = transf(amenities) * cf_wage / cf_c_price / 100 * 1e4
  ) %>%
  group_by(nonattainment) %>%
  summarise(
    total_income_gain = sum(cf_labor * total_wage_gain * workers_total, na.rm = T),
    cons_income_gain = sum(cf_labor * cons_wage_gain * workers_total, na.rm = T),
    amenity_income_gain = sum(cf_labor * amenity_wage_gain * workers_total, na.rm = T),
    wages = sum(cf_wage / cf_c_price * 1e4 * workers_total * cf_labor, na.rm = T)
  ) |>
  mutate(total_income_gain = cons_income_gain + amenity_income_gain)

welfare_base_first <- fread("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  left_join(nonattainment) %>%
  select(fips, welfare, cf_labor, cf_wage, real_wage, amenities, cf_c_price, nonattainment) %>%
  mutate(
    total_wage_gain = transf(welfare) * cf_wage / cf_c_price / 100 * 1e4,
    cons_wage_gain = transf(real_wage) * cf_wage / cf_c_price / 100 * 1e4,
    amenity_wage_gain = transf(amenities) * cf_wage / cf_c_price / 100 * 1e4
  ) %>%
  group_by(nonattainment) %>%
  summarise(
    total_income_gain = sum(cf_labor * total_wage_gain * workers_total, na.rm = T),
    cons_income_gain = sum(cf_labor * cons_wage_gain * workers_total, na.rm = T),
    amenity_income_gain = sum(cf_labor * amenity_wage_gain * workers_total, na.rm = T),
    wages = sum(cf_wage / cf_c_price * 1e4 * workers_total * cf_labor, na.rm = T)
  ) |>
  mutate(total_income_gain = cons_income_gain + amenity_income_gain)

welfare_base_pe <- fread("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-off_mobility-off_mktacc-off_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  left_join(nonattainment) %>%
  select(fips, welfare, cf_labor, cf_wage, real_wage, amenities, cf_c_price, nonattainment) %>%
  mutate(
    total_wage_gain = transf(welfare) * cf_wage / cf_c_price / 100 * 1e4,
    cons_wage_gain = transf(real_wage) * cf_wage / cf_c_price / 100 * 1e4,
    amenity_wage_gain = transf(amenities) * cf_wage / cf_c_price / 100 * 1e4
  ) %>%
  group_by(nonattainment) %>%
  summarise(
    total_income_gain = sum(cf_labor * total_wage_gain * workers_total, na.rm = T),
    cons_income_gain = sum(cf_labor * cons_wage_gain * workers_total, na.rm = T),
    amenity_income_gain = sum(cf_labor * amenity_wage_gain * workers_total, na.rm = T),
    wages = sum(cf_wage / cf_c_price * 1e4 * workers_total * cf_labor, na.rm = T)
  ) |>
  mutate(total_income_gain = cons_income_gain + amenity_income_gain)

# first best is relative to 1997 nonattainment so need to difference against 1997 nonattainment relative to no nonattainment
welfare_base_first[2:4] <- (welfare_base_nonattain)[2:4] - welfare_base_first[2:4]
welfare_first_industry[2:4] <- (welfare_base_industry)[2:4] - welfare_first_industry[2:4]

base_welfare_table <- tribble(
  ~type, ~total_percent, ~total_dollars, ~amenities_percent, ~amenities_dollars, ~prod_percent, ~prod_dollars,
    "Aggregate", welfare_base$total, sum(welfare_base_nonattain$total_income_gain), welfare_base$amenity_total, sum(welfare_base_nonattain$amenity_income_gain), welfare_base$consumption_total, sum(welfare_base_nonattain$cons_income_gain),
  "Manufacturing", welfare_base$manufacturing, welfare_base_industry$total_income_gain[1], welfare_base$manu_amenity, welfare_base_industry$amenity_income_gain[1], welfare_base$manu_cons, welfare_base_industry$cons_income_gain[1],
  "Nonmanufacturing", welfare_base$other, welfare_base_industry$total_income_gain[2], welfare_base$other_amenity, welfare_base_industry$amenity_income_gain[2], welfare_base$other_cons, welfare_base_industry$cons_income_gain[2],
  "Nonemployed", welfare_base$unemp, NA, welfare_base$unemp_amenity, NA, NA, NA,
  "Attainment Counties", welfare_base$welfare_attain, welfare_base_nonattain$total_income_gain[1], welfare_base$amenity_attain, welfare_base_nonattain$amenity_income_gain[1], welfare_base$cons_attain, welfare_base_nonattain$cons_income_gain[1],
  "Nonattainment Counties", welfare_base$welfare_nonattain, welfare_base_nonattain$total_income_gain[2], welfare_base$amenity_nonattain, welfare_base_nonattain$amenity_income_gain[2], welfare_base$cons_nonattain, welfare_base_nonattain$cons_income_gain[2],
  "Aggregate", welfare_pe$total, sum(welfare_base_pe$total_income_gain), welfare_pe$amenity_total, sum(welfare_base_pe$amenity_income_gain), welfare_pe$consumption_total, sum(welfare_base_pe$cons_income_gain),
  "Manufacturing", welfare_pe$manufacturing, welfare_pe_industry$total_income_gain[1], welfare_pe$manu_amenity, welfare_pe_industry$amenity_income_gain[1], welfare_pe$manu_cons, welfare_pe_industry$cons_income_gain[1],
  "Nonmanufacturing", welfare_pe$other, welfare_pe_industry$total_income_gain[2], welfare_pe$other_amenity, welfare_pe_industry$amenity_income_gain[2], welfare_pe$other_cons, welfare_pe_industry$cons_income_gain[2],
  "Nonemployed", welfare_pe$unemp, NA, welfare_pe$unemp_amenity, NA, NA, NA,
  "Attainment Counties", welfare_pe$welfare_attain, welfare_base_pe$total_income_gain[1], welfare_pe$amenity_attain, welfare_base_pe$amenity_income_gain[1], welfare_pe$cons_attain, welfare_base_pe$cons_income_gain[1],
  "Nonattainment Counties", welfare_pe$welfare_nonattain, welfare_base_pe$total_income_gain[2], welfare_pe$amenity_nonattain, welfare_base_pe$amenity_income_gain[2], welfare_pe$cons_nonattain, welfare_base_pe$cons_income_gain[2],
  "Aggregate", welfare_first$total, sum(welfare_base_first$total_income_gain), welfare_first$amenity_total, sum(welfare_base_first$amenity_income_gain), welfare_first$consumption_total, sum(welfare_base_first$cons_income_gain),
  "Manufacturing", welfare_first$manufacturing, welfare_first_industry$total_income_gain[1], welfare_first$manu_amenity, welfare_first_industry$amenity_income_gain[1], welfare_first$manu_cons, welfare_first_industry$cons_income_gain[1],
  "Nonmanufacturing", welfare_first$other, welfare_first_industry$total_income_gain[2], welfare_first$other_amenity, welfare_first_industry$amenity_income_gain[2], welfare_first$other_cons, welfare_first_industry$cons_income_gain[2],
  "Nonemployed", welfare_first$unemp, NA, welfare_first$unemp_amenity, NA, NA, NA,
    "Attainment Counties", welfare_first$welfare_attain, welfare_base_first$total_income_gain[1], welfare_first$amenity_attain, welfare_base_first$amenity_income_gain[1], welfare_first$cons_attain, welfare_base_first$cons_income_gain[1],
  "Nonattainment Counties", welfare_first$welfare_nonattain, welfare_base_first$total_income_gain[2], welfare_first$amenity_nonattain, welfare_base_first$amenity_income_gain[2], welfare_first$cons_nonattain, welfare_base_first$cons_income_gain[2],

)

base_welfare_table <- base_welfare_table %>%
  mutate(
    total_dollars = total_dollars,
    amenities_dollars = amenities_dollars,
    prod_dollars = prod_dollars
  ) %>%
  mutate(
    across(ends_with("percent"), ~ round(.x, 2)),
    across(ends_with("dollars"), ~ round(.x, 0)),
    across(total_percent:prod_dollars, ~ as.character(.x)),
    across(total_percent:prod_dollars, ~ ifelse(is.na(.x), "-", .x))
  )

base_table_out <- kbl(base_welfare_table,
  format = "latex",
  caption = "Welfare impacts of nonattainment in 1997.",
  label = "base_welfare",
  booktabs = T,
  digits = 1,
  escape = FALSE,
  linesep = c(""),
  col.names = linebreak(c("", "\\%", "Billion \\$", "\\%", "Billion \\$", "\\%", "Billion \\$"), align = "c"),
  align = c("l", rep("c", dim(base_welfare_table)[2] - 1))
) %>%
  add_header_above(c(" " = 1, "Total" = 2, "Amenity" = 2, "Consumption" = 2)) %>%
  kable_styling(latex_options = "scale_down") %>%
  pack_rows(index = c("A. 1997 Nonattainment" = 6, "B. No Economic/Physical Geography" = 6, "C. First-Best Emissions Pricing" = 6)) %>%
  # row_spec(c(1), italic = T) %>%
  footnote(
    general = "Welfare is computed as the equivalent variation of (A) the observed nonattainment status in 1997, (B) the observed nonattainment status in 1997 if there is no labor reallocation, market access is held fixed, and pollution damages are only within the emitting county, or (C) first-best emissions pricing, relative to a counterfactual in which no counties are in nonattainment or face emissions pricing. The simulations in (A) and (C) account for labor reallocation, trade, and atmospheric transport of pollution. The first-best result sets the optimal location-specific nonnegative emission prices. Numbers may not sum up fully due to rounding.",
    threeparttable = T,
    fixed_small_size = T
  )

fileConn <- file("output/table_3_welfare_table.tex")
writeLines(base_table_out, fileConn)
close(fileConn)

##########################
##### ROBUSTNESS TABLE
##########################

welfare_base <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "Base Model")
welfare_iota_l <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.274_gamma-0.481_iota-0.495_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\iota = 0.495$")
welfare_iota_h <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.274_gamma-0.481_iota-2.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\iota = 2$")
welfare_theta_h <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-8_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\theta = 8$")
welfare_theta_l <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-2.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\theta = 2$")
welfare_alpha_h <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.548_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\alpha = 0.548$")
welfare_alpha_l <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.137_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\alpha = 0.137$")
welfare_gamma_h <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.274_gamma-1_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\gamma = 1$")
welfare_gamma_l <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.274_gamma-0.2405_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\gamma = 0.2405$")
welfare_xi_l <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-0.5_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\xi$s Halved")
welfare_xi_h <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-2_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\xi$s Doubled")
welfare_xi_hh <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-4_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "$\\xi$s Quadrupled")
welfare_cong <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-on_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "Add Congestion/Agglomeration")
welfare_vsl <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-inc.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "VSL Increases in Income")
welfare_p3 <- fread("data/counterfactual/welfare_total_opt-base_prod-0.03_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "Increase Productivity 3 Percent")
welfare_m3 <- fread("data/counterfactual/welfare_total_opt-base_prod--0.03_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:other_amenity, ~ transf(.x))) %>%
  mutate(type = "Decrease Productivity 3 Percent")

welfare_table <- bind_rows(
  welfare_base, welfare_theta_l, welfare_theta_h, welfare_alpha_l, welfare_alpha_h, welfare_gamma_l, welfare_gamma_h, welfare_iota_l, welfare_iota_h, welfare_xi_l, welfare_xi_h, welfare_xi_hh, welfare_cong, welfare_vsl, welfare_p3, welfare_m3
) %>%
  select(type, total, manufacturing, other, welfare_attain, welfare_nonattain) %>%
  mutate(
    across(2:6, ~ round(.x, 2)),
    across(2:6, ~ as.character(.x))
  )

base_table_out <- kbl(welfare_table,
  format = "latex",
  caption = "Welfare impacts of nonattainment in 1997 under different model parameters or structural assumptions.",
  label = "base_welfare_robust",
  booktabs = T,
  digits = 2,
  linesep = c(""),
  escape = FALSE,
  col.names = linebreak(c("Change", "Total \\%", "Manuf. \\%", "Nonmanuf. \\%", "Attain. \\%", "Nonattain. \\%"), align = "c"),
  align = c("lccccc")
) %>%
  kable_styling() %>%
  # row_spec(c(1), bold = T) %>%
  pack_rows(index = c("No Changes" = 1, "A. Parameter Changes" = 11, "B. Structural Changes" = 4)) %>%
  footnote(
    general = "Welfare is computed as the equivalent variation of the observed nonattainment statuses in 1997 relative to a counterfactual where no counties are in nonattainment. Each row only makes one change relative to the base model in the first row.",
    threeparttable = T,
    fixed_small_size = T
  )

fileConn <- file("output/table_d2_welfare_table_robust.tex")
writeLines(base_table_out, fileConn)
close(fileConn)
